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Abstract. We study a general class of infimal convolution type regularisation functionals 
suitable for applications in image processing. These functionals incorporate a combination 
of the total variation (TV) seminorm and norms. A unified well-posedness analysis is 
presented and a detailed study of the one dimensional model is performed, by computing 
exact solutions for the corresponding denoising problem and the case p = 2. Furthermore, 
the dependency of the regularisation properties of this infimal convolution approach to the 
choice of p is studied. It turns out that in the case p = 2 this regulariser is equivalent 
to Huber-type variant of total variation regularisation. We provide numerical examples for 
image decomposition as well as for image denoising. We show that our model is capable of 
eliminating the staircasing effect, a well-known disadvantage of total variation regularisation. 
Moreover as p increases we obtain almost piecewise affine reconstructions, leading also to a 
better preservation of hat-like structures. 

Keywords: Total Variation, Infimal convolution, Denoising, Staircasing, norms. Image 
decomposition 


1. Introduction 

In this paper we introduce a family of novel TV-L^ infimal convolution type functionals 
with applications in image processing: 

(1.1) Ju) inf a\\Du-w\\M +l3\\w\\hP(Q)^ Q!,/3>0 and p>l. 

Here || • \\m denotes the Radon norm of a measure. The functional (1.1) is suitable to be used as 
a regulariser in the context of variational non-smooth regularisation in imaging applications. 
We study the properties of (1.1), its regularising mechanism for different values of p and apply 
it successfully to image denoising. 

1.1. Context. After the introduction of the total variation (TV) for image reconstruction 
purposes [27], the use of non-smooth regularisers has become increasingly popular during the 
last decades (cf. [8]). They are typically used in the context of variational regularisation, 
where the reconstructed image is obtained as a solution of a minimisation problem of the 
type: 

(1.2) + ^(m). 

The regulariser is denoted here by T. We assume that the data /, defined on a domain 
C M^, have been corrupted through a bounded, linear operator T and additive (random) 
noise. Different values of s can be considered for the first term of (1.2), the fidelity term. 
For example, models incorporating a fidelity term (resp. L^) have be shown to be efficient 
for the restoration of images corrupted by Gaussian noise (resp. impulse noise). Of course, 
other types of noise can also be considered and in those cases the form of the fidelity term 
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is adjusted accordingly. Typically, one or more parameters within T balance the strength of 
regularisation against the fidelity term in the minimisation (1.2). 

The advantage of using non-smooth regularisers is that the regularised images have sharp 
edges (discontinuities). For instance, it is a well-known fact that TV regularisation promotes 
piecewise constant reconstructions, thus preserving discontinuities. However, this also leads 
to blocky-like artifacts in the reconstructed image, an effect known as staircasing. Recall at 
this point that for two dimensional images u G L^(fl), the definition of the total variation 
functional reads: 


TV(m) sup<^ / udiYcfxix : v € ||0||^ < i 


(1.3) 

The total variation uses only first-order derivative information in the regularisation process. 
This can be seen from that fact that for TV {u) < oo the distributional derivative Du is a finite 
Radon measure and TV(?i) = Moreover if G W^’^(fl) then TV(?i) = |Vi/| dx, 

i.e., the total variation is the norm of the gradient of u. Higher-order extensions of the 
total variation functional are widely explored in the literature e.g. [4,5,9,11,12,20-22,25]. The 
incorporation of second-order derivatives is shown to reduce or even eliminate the staircasing 
effect. The most successful regulariser of this kind is the second order total generalised 
variation (TGV) introduced by Bredies et al. [5]. Its definition reads 


(1.4) 


TGV^ o(u) := min a\\Du — w\\m -\- P\\£w\\m- 
weBB{n) 


Here a, [3 are positive parameters and BD(fl) is the space of functions of bounded deformation, 
i.e., the space of all L^(fl) functions u;, whose symmetrised distributional derivative £w is a 
finite Radon measure. This is a less regular space than the usual space of functions of bounded 
variation BV(fl) for which the full gradient Du is required to be a finite Radon measure. Note 
that if the variable w in the definition (1.4) is forced to be the gradient of another function 
then we obtain the classical infimal convolution regulariser of Chambolle-Lions [9]. In that 
sense TGV can be seen as a particular instance of infimal convolution, optimally balancing 
first and second-order information. 

In the discrete formulation of TGV (as well as for TV) the Radon norm is interpreted as 
an norm. The motivation for the current and the follow-up paper is to explore the capa¬ 
bilities of norms within first-order regularisation functionals designed for image processing 
purposes. The use of norms for p > 1 has been exploited in different contexts - infinity 
and p-Laplacian (cf. e.g. [14] and [19] respectively). 

1.2. Our contribution. Comparing the definition (1.1) with the definition of TGV in 
(1.4), we see that the Radon norm of the symmetrised gradient of w has been substituted by 
the norm of u;, thus reducing the order of regularisation. Up to our knowledge, this is the 
first paper that provides a thorough analysis of TV-L^ infimal convolution models (1.1) in this 
generality. We show that the minimisation in (1.1) is well-defined and that T\L^ ^{u) < oo 
if and only if TV(i/) < cx). Hence TVL^^ regularised images belong to BV(fl) as desired. 

In order to get more insight in the regularising mechanism of the TVL^ ^ functional we 
provide a detailed and rigorous analysis of its one dimensional version of the corresponding 
fidelity denoising problem 
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^ll/-^llL2(a) + TVLly«). 




(1.5) min , 

nGBV(O) 2 

For the denoising problem (1.5) with p = 2 we also compute exact solutions for simple one 
dimensional data. We show that the obtained solutions are piecewise smooth, in contrast to 
TV (piecewise constant) and TGV (piecewise affine) solutions. Moreover, we show that for 
p = 2, the 2-homogeneous analogue of the functional (1.1) 


(1.6) 


F(u) = min a\ 

i(;GL2(0) 


Du 


w 


\m + ^I|w^IIl2(q)) 
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is equivalent to a variant of Huber TV [17], with the functional (1.6) having a close connection 
with (1.1) itself. Huber total variation is a smooth approximation of total variation and even 
though it has been widely used in the imaging and inverse problems community, it has not 
been analysed adequately. Hence, as a by-product of our analysis, we compute exact solutions 
of the one dimensional Huber TV denoising problem. 

We proceed with exhaustive numerical experiments focusing on (1.5). Our analysis is con¬ 
firmed by the fact that the analytical results coincide with the numerical ones. Furthermore, 
we observe that even though a first-order regularisation functional is used, we are capable of 
eliminating the staircasing effect, similarly to Huber TV. By Bregmanising our method [23], 
we are also able to enhance the contrast of the reconstructed images, obtaining results very 
similar in quality to the TGV ones. We observe numerically that high values of p promote 
almost affine structures similar to second-order regularisation methods. We shed more light 
of this behaviour in the follow-up paper [7] where we study in depth the case p = oc. Let 
us finally note that we also consider a modified version of the functional (1.1) where w is re¬ 
stricted to be a gradient of another function leading to the more classical infimal convolution 
setting. Even though, this modified model is not so successful in staircasing reduction, it is 
effective in decomposing an image into piecewise constant and smooth parts. 

1.3. Organisation of the paper. After the introduction we proceed with the introduction 
of our model in Section 2. We prove the well-posedness of (1.1), we provide an equivalent 
definition and we prove its Lipschitz equivalence with the TV seminorm. We finish this 
section with a well-posedness result of the corresponding TVL^ ^ regularisation problem using 
standard tools. 

In Section 3 we establish a link between the TVL^^ functional and its p-homogeneous 
analogue (using the p-th power of || • ||lp(0))- The p-homogeneous functional (for p = 2) is 
further shown to be equivalent to Huber total variation. 

We study the corresponding one dimensional model in Section 4 focusing on the fidelity 
denoising case. More specifically, after deriving the optimality conditions using Fenchel- 
Rockafellar duality in Section 4.1, we explore the structure of solutions in Section 4.2. In 
Section 4.3 we compute exact solutions for the case p = 2, considering a simple step function 
as data. 

In Section 5 we present a variant of our model suitable for image decomposition purposes, 
i.e., geometric decomposition into piecewise constant and smooth structures. 

Section 6 focuses on numerical experiments. Confirmation of the obtained one dimensional 
analytical results is done in Section 6.2, while two dimensional denoising experiments are 
performed in Section 6.3 using the split Bregman method. There, we show that our approach 
can lead to elimination of the staircasing effect and we also show that by using a Bregmanised 
version we can also enhance the contrast, achieving results very close to TGV, a method 
considered state of the art in the context of variational regularisation. We finish the section 
with some image decomposition examples and we summarise our results in Section 7. 

In the appendix, we remind the reader of some basic facts from the theory of Radon 
measures and BV functions. 

2. Total variation and regularisation 

In this section we introduce the TV-L^ functional (1.1) as well as some of its main prop¬ 
erties. For (T, /3 > 0 and 1 < p < oo, we define TVL^ ^ : L^(fl) ^ M as follows: 

TVL'„(„) 

- w\\m + IMhPiii) ■ 

The next proposition asserts that the minimisation in (1.1) is indeed well-defined. We omit 
the proof, which is based on standard coercivity and weak lower semicontinuity techniques: 

Proposition 2.1. Let u G BV(fl) with 1 < p < oc and d,/? > 0. Then the minimum in the 
definition (1.1) is attained. 
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Another useful formulation of the definition (1.1) is the dual formulation: 

(2.1) TVLl^iu) = supS^j^udw(pdx : (p e Cl{n),\\(p\\^ < a, ||^||L.(n) < ^| , 

The next Proposition shows that the two expressions coincide indeed. 

Proposition 2.2. Let u G BV(fl) and 1 < p < oo then 

in a||Il«-Ty||_y^+^||w;||LP(m =sup<^ / udwcpdx : 0 G ||^||^ < a, \\(p\\i^,m </3 

Un 


min 

W^LP{Q) 


Proof. First notice that in (2.1), we can replace C^(fl) by CQ(fl), since Cl{Pt) = ^^(fl) with 
the closure taken with respect to the uniform norm. We define 

Fi:X with Fiicp) = I{||.||L,(f,)</3}(<?^)’ 

F2 : X ^ R, with F2{4>) = I|||.|| — J udiv(pdx. 

Then, we can rewrite (2.1) as 


TVL',/«) = - 


inf 

(f)ex 


I — I u dive/) dx \ = — inf Fi (0) + F 2 (0) • 

[ Jq J 


W^Wlq (n)—h 

The Fenchel-Rockafellar duality theory, see [13], allows to establish a relation between the 
primal problem 


and its dual 


■ inf Fi(0) + F 2 (</.), 

0GA 


min F^{-w) + F 2 iw). 

w^X* 


Here Ff and denote the convex conjugate of Fi and F 2 respectively. In order to obtain 
such a connection, we follow [2] where it suffices to show that 

X (domF 2 — domFi) 

A>0 

is a closed vector space. Indeed, we have that 

X (domF 2 — domFi) C X, 

A>0 

and for every 0 G X, we can write 0 = X{X~^(p — 0) with ^ cind 0 G domFi. 

Hence, IJ A (domF 2 — domFi) = X is a closed vector space and there is no duality gap i.e., 
A>0 

inf {Fi(0) + F^m + min {Ff{w) + F^{w)} = 0. 

0GA k;GA* 


Finally, we have 


and similarly. 


FU-w) = sup {(l),w) = 13 ||w||lp(q) , 

4>ec^(n) 

F 2 {w) = sup {w^ (j)) + {u^ 0^) = a \\Du — w\\j^ . 


Thus the desired equality is proven. 


□ 
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Remark 2.3. The dual formulation o/TVL^^ : M. is useful sinee one ean easily 

derive that TVL^ ^ is lower semieontinuous with respeet to the strong topology sinee it is 
a pointwise supremum of eontinuous funetions. 

The following lemma shows that the TVL^ ^ functional is Lipschitz equivalent to the total 
variation seminorm. 

Lemma 2.4. Let u G and 1 < p < oo. Then TVL^ ^(?i) < oc if and only if u ^ BV(fi) 

and there exist eonstants 0 < (7i,(72 < oc sueh that 

(2.2) C 2 \\Du\\j^ < TVL^,^(«) < Cl \\Du\\j^ . 

Proof. Let u G BV(fi). Using (1.1) we have that 

TVL^,/3(«) < a \\Du - w\\_^ + p . 

for every w G Setting w = 0 and Ci = a, we obtain 

TVL^,yu) < Cl \\Du\\j^ . 

For the other direction, we have for any w € L^’(f2) C L^(f2) by the triangle inequality 

1 

\\Du\\^ < \\Du - w\\^ + ||w||Li(n) < \\Du - w\\^ + 

< C(||F>n-n;||_y^( + ||«^||LP(f2))> 

1 

with C = max(l, For appropriate choice of C 2 we obtain 

C 2 \\Du\\j^ < Q;||F>n '^Wm II'^IIlp(O) ’ 

which yields the left-hand side inequality. □ 


Having shown the basic properties of the TVL^ ^ functional, we can use it as a regulariser 
for variational imaging problems, by minimising 


(2.3) 


1 

mm - 
nGBV(O) S 


\\f-Tu\\ls^^^+TYLl^{u)^ 


s > 1, 


where T : L^(fl) ^ L^(fl) is a bounded, linear operator and / G L^(fl). We conclude our 
analysis with existence and uniqueness results for the minimisation problem (2.3). 


Theorem 2.5. Let 1 < p < 00 and f G L^(fl). //T(Tb) ^ 0 then there exists a solution 
u G L^(fl) n BV(fl) for the problem (2.3). If s > 1 and T is injeetive then the solution is 
unique. 


Proof. The proof is a straightforward application of the direct method of calculus of variations. 
We simply take advantage of (2.2) and the compactness theorem in BV(fl) along with the 
lower semicontinuity property of TVL^ We also refer the reader to the corresponding proofs 
in [25,28]. ’ □ 


Since we are mainly interested in studying the regularising properties of TVL^ from now 
on we focus on the case where 5 = 2 and T is the identity function (denoising task) where 
rigorous analysis can be carried out. We thus define the following problem 

or equivalently 

1 2 

(V) min - 11/ - n||L2(m + aWDu - w\\m + P\\w\\LP{n)- 

nGBV(S2j Z ^ ^ 

weLP{n) 
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3. The ^-homogeneous analogue and relation to Huber TV 

Before we proceed to a detailed analysis of the one dimensional version of (P), in this 
section we consider its p-homogeneous analogue 

(Vp-hom) -u\\l2(^^) + a\\Du-w\\M +I < p < oo. 

weLP{Q) 

We show in Proposition 3.2 that there is a strong connection between the models (V) and 
(Vp-hom)- The reason for the introduction of Vp-hom is that, in certain cases, it is technically 
easier to derive exact solutions for {Vp-hom) rather than for {V) straightforwardly, see Section 
4.3. Moreover, here we can guarantee the uniqueness of the optimal since 

= argmin a \\Du — w\\j^ + — H'R'IlLPfo) ? 
weLP{Q) P ^ ^ 

and thus is unique as a minimiser of a strictly convex functional since 1 < p < oo. Hence, 
compared to ('P), an optimal solution pair of {Vp-hom) is unique. The next proposition says 
that, unless / is a constant function then the optimal w in {Vp-hom) cannot be zero but 
nonetheless converges to 0 as /3 ^ oc. In essence, this means that one cannot obtain TV type 
solutions with the p-homogeneous model. 

Proposition 3.1. Let 1 < p < oo, / G and let be an optimal solution pair of 

the p-homogeneous problem {Vp-hom)- Then rt;* = 0 i/ and only if f is a eonstant funetion. 
For general data /, we have that ^ 0 in for [3 oo. 


Proof. It follows immediately that if / is constant then (0, /) is the optimal pair for {Vp-hom)- 
Suppose now that (n;*,n*) solve {Vp-hom)- If is easy to check that the following also hold: 

(3.1) w* e argmin a\\Vu* - w||Lim) + -\\w\\ip(n)^ 

weLP(Q) P ^ ’ 

(3.2) = argmin Lf - «||l 2 (q) + a\\Du - w*\\m- 

«gbv(Q) ^ 

In particular, (3.1) implies that 

(3.3) 0 G ad\\Vu* - •||Li(f 2 )(«^*) + P\w*\P-^w*, 

Suppose now that n;* = 0. Then (3.2) becomes 

(3.4) u* ^ argmin h\f - u\\l 2 rm + a\\Du\\M- 

nGBV(O) ^ 

Furthermore, since (0,n*) solve {Vp-hom)^ then for every u G C^{fl) and e > 0, the pair 
(eVn,n* -h eu) is suboptimal for {Vp-hom), i-c., 

^11/ - «*IIl 2 (Q) + a\\Du*\\M < i||(/ - U*) - €«||^ 2 (q) + a\\D{u* + eu) - eVu\\M + 
from which we take 

^11/ - ’^*llL2(n) < \\\if - u*) - e’^llL2(n) + ^lle'^^IIL(Q) 

0 < j\\u\\l2^a) ~ u*)udx + ^hSIu\\l^^^y 

By dividing the last inequality by e and taking the limit e ^ 0 we have that j^{f — u')u dx < 0. 
By considering the analogous perturbations n* — en , we obtain similarly that J^{f — u'')u dx > 
0 and thus 

/ (/ - u^)udx = 0, Vn G C^{n). 

Jq 
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Hence u" — f and since u" solves (3.4) this can only happen when / is a constant function. 

For the last part of the proposition, (supposing / ^ 0), simply observe that for every 
u G BV(i4) and w G L^(i4) we have that 


1 

2 


11/-^* 


||2 

liL2(a) 


+ a\\Du* 




'^ll At + 


p 


w 


IIP 

IIlp(O)’ 


and setting = ip = 0, we obtain 


1 

P 


\w 


IIP ^ J_ 

IIlp(O) - 2/3 


II/IIl2(o)’ 


and thus ^ 0 when /3 ^ oo. 


□ 


We can further establish a connection between the 1-homogeneous (V) and the p-homogeneous 
model {'Pp—hom)' 

Proposition 3.2. Let 1 < p < oo and f G L^(ri) not a constant. A pair (ip*,?i*) is a 
solution of (Vp-hom) parameters {a^ {dp^hom) if only if it is also a solution of (V) 
with parameters {a.Pi^hom) where Pi-hom = Pp-hom\\w''\\^^^^^y 


Proof Since / is not a constant by the previous proposition we have that rp* ^ 0. Note that 
for an arbitrary function u G BV(i4): 


w* G argmin a \\Du — w\\ ^ 

weLPiQ) P ^ 

0 ead\\Du - -11^ (ip*) + 

fdi- 


0 ^ ad \\Du — 


lAt 


(^*) + 




hom I *|»—2 * ^ 

-— \w p W 

||P-1 I I 

IIlp(O) 


w* € argmin a \\Du - w\\_f^ + Pi-hom II«^IIlp(Q) • 

weLP{Q) 


This means that ip* is an admissible solution for both problems (V) and {Vp-hom)^ with the 
corresponding set of parameters {a^ f3i-hom) si^nd {a, f3p-hom) respectively. The fact that the 
same holds for u* as well, comes from the fact that in both problems we have 

u* e argmini||/-n||L 2 m) + a||-Dn-ru*||x. 

«gbv(Q)^ 


□ 


Finally, it turns out that for p = 2, problem {Vp-hom) is essentially equivalent to the widely 
used Huber total variation regularisation, [17]. We show that in the next proposition. 


Proposition 3.3. Consider the functional F : BV(fl) ^ M with 
(3.5) Fiu)= ^}^^a\\Du - w\\m + ^\\w\\l2ra). 


Then 


where p : ^ M with 


F{u)= [ p{Vu)dxFa\D^u\{n), 

Jn 


(p{x) 


{a\x\-^, kl > f) 

I ^ I ^ 12 I I — 

f 2 
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Proof. We have 


F{u) = min a\\Du-w\\M + ^\Mh(Q) 


(3 


= a\D^u\(Q)-\-a min / \Vu — w\-\-^\w\‘^ dx. 

weL‘^(fi)jQ ‘2a 


So we focus on the minimisation problem 


(3.6) 


/3 


min / iVu — w\-\-—\wr dx. 
weL^{Q)jQ 2a 


Baring in mind that (as it can easily checked) for c G and A > 0, 

tie 


A, 


if |c| > 


argmin |c-y| +-|yp = <^ ^1"=! 

^ c if|c|<^, 


and 




it is straightforwardly verified setting \ — ^ja that the function 

^ a Vu 


w — 


belongs to C and solves (3.6) with optimal value equal to A (p(Vu) dx. 


>« \ + \/uXf 


□ 


4. The one dimensional case 

In order to get more insight into the structure of solutions of the problem (P), in this 
section we study its one dimensional version. As above, we focus on the finite p case, i.e., 
1 < p < oc. The case p = oo leads to several additional complications and will be subject 
of a forthcoming paper [7]. For this section 14 C M is an open and bounded interval, i.e., 
14 = (a, 6). Our analysis follows closely the ones in [6] and [24] where the one dimensional 
L^-TGV and L^-TGV problems are studied respectively. 

4.1. Optimality conditions. In this section, we derive the optimality conditions for the 
one dimensional problem (V). We initially start our analysis by defining the predual problem 
(P*), proving existence and uniqueness for its solutions. We will employ again the Fenchel- 
Rockafellar duality theory in order to find a connection between their solution pairs. 

We define the predual problem (P*) as {q Holder conjugate to p): 

{V*) - inf fc^'dx + i l^{<f>'?dx : <f> € HS(fi), 11011^,(0) < P, U\L < «} • 

Existence and uniqueness can be verified by standard arguments: 

Proposition 4.1. For f G L?{Pl), the predual problem (P*) admits a unique solution in 

Observe now that we can also write down the predual problem (P*) using the following 
equivalent formulation: 


- inf Fi( 0 ,O + F2(i^(0,e), 

(<l>,n)ex 
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(4.1) 



where X = Hj(fi) x y = Hj(fi) x L2(fi) and 

K:X^Y, i^(</.,O = (^-0,O, 

(4.2) F,:X^W, with Fi(0,O = I{|M|,,(,^</3}(<^) + I{|MU<o}(0, 

F 2 :Y ^R, with F 2 (</>,'*/’) = I{o}(</’) + [ f'’P dx + ^ f ip'^dx. 

Jn ^ Jn 

We denote the infimum in (P*) as infP*. Then, it is immediate that 

- inf iP* = - inf Fi{(I>,0 + F 2 {K{(I>, 0 ). 

The dual problem of (4.1), see [13], is defined as 


(4.3) 


min Fi{—K'^{w^ u)) + F 2 {w^ u)^ 
{w,u)eY* 


where here denotes the adjoint of K. Let (cr, r) be elements of HQ(ri)* x HQ(fi)* acting as 
distributions. For the convex conjugate of Fi, we write 

(4.4) Fi((T,r)= sup (cr, 0) + (r, 0 ^ sup {a,(f))+a sup (t,^ • 

?gh1(Q) 

ll^llM(n)</3 ll0llL9(n)<l IISllco<l 

ll€lloo<« 

However, by standard density arguments we have: 

(4.5) Fi(<T,r) = ^ sup {a,(l))+a sup (r, ^ ^ IkllLPm) +«IklU • 

ll^llL9(fi)—^ IlCiloo —^ 

Moreover, let : T* ^ X* with 

{-K*{w, u), (</., 0) = - {{w, u), Ki<f>, 0) =-{{w,u),{^- cf>, e')> 

= - {w, 0 + {w, <t>) - {u, O = (dDu -w,^) + {w, (f)) . 

Hence, we obtain 

(4.6) Fa-K*iw, u)) = P ||u;||L.(n) + a\\Du- w\\j^ , 
and 

F 2 {w,u) = sup {w,(j)) + {u,'lp) - (/,V’) - X [ V’^ 

i4>,^)&Y ^ Jn 

A=o 


1 


= sup {u- f, Ip) - - / pj^dx 

(4.7) ^=(i|MlhQ,)*(«-/) = i|l«-/llhQ,. 

Therefore, we have proved the following: 

Proposition 4.2. The dual problem of (F*) is equivalent to the primal problem (V) in the 
sense that (w^u) G T* is a solution of the dual of (F*) if and only if (w^u) G x BV(ri) 

is a solution of (V). 


It remains to verify that we have no duality gap between the two minimisation problems 
{V) and (F*). The proof of the following proposition follows the proof of the corresponding 
proposition in [6]. We slightly modify for our case. 

Proposition 4.3. Let Fi^F 2 ^K be defined as in (4.2). Then 

(4.8) ^ = U ^(domFa - i^(domFi)) 

A>0 
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is a closed vector space and thus [2] 

(4.9) mm F,{{4>,0) + F2{K{4>,0) + ^ min Fr-K*{w,u)) + F^{{w,u)) = 0. 

(w,u)eY* 

Proof. Let (0, and define f^o{x) = C 1 +C 2 X, where ci, C 2 are constants that are uniquely 

determined by the following conditions 

(4.10) / dx = / fj{x)dx^ / xfjo{x) dx = / xfj{x) + (/){x) dx. 

Jn Jn Jn Jn 

Let ^(x) = J^('0o~'0)(2/) dy> Since by construction, = '0o~'0 ^ L^(fl) with ^{a) — ^{b) — 0, 
we have that ^ G HQ(fl). Furthermore, let (j) = —(0 + C) ^ and (0,C) ^ ^ with 

(^^V’) = (C - ^,V’o - ^') = (0,V’o) - - 0,0 = (O,0o) - K{(f),Cj- 

Choosing appropriately A > 0 such that IIA 


- 1 , 


lm(o) 


< /3, A < (a, we can write 


(0,0) = A((0, A-Vo) - i^(A-V, A-'C), 

with domF 2 = {0} x L2(fi) and domFi = {(0,0 : II0IIl9(Q) ^ IIOL ^ V- Since (0,0) € 
were chosen arbitrarily, (4.8) holds. 

Since there is no duality gap, we can find a relationship between the solutions of (P*) and 
{V) via the optimality conditions, see [13, Prop. 4.1 (III)]. 

Theorem 4.4 (Optimality conditions). Let 1 < p < oo and f G L^(fl). A pair (w^u) G 
L^(54) X BV(fl) is a solution of (V) if and only if there exists a function (j) G HQ(fl) such that 

0' = ^ - /, 

(j) G aSgn{Du — w)^ 


(4.11) 

and 


(4.12) 


0 G |0 G L9(f^) : 


L'3(a) 


</3 


ifw = 0, 
if w ^ 0. 


Ilp(q) >a\\Du-w\\j^ + (3 ||u;||lp(q) + {{a - n;), 0) + (r - {Du -w),^) ^ 


Proof. Since there is not duality gap, the optimality conditions read [13, Prop. 4.1 (III)]: 

(4.13) {(l>,0edFf{-K*{w,u)), 

(4.14) K{f>,0^dF^{w,u), 

for every (0,0 {w,u) solutions of (T*) and (P) respectively. Note that in dimension one 

we have HQ(n) C (70(12). Hence, for every (a, r) G X*, we have the following: 

Fr(T,T) > Fr-K*{w,u)) + {{a,T) + K*{w,u),{cp,0) ^ 

« VI V + 0 IVIlLP(a) ^ - M\m + 0 I VllLP(a) + (V 0 - i-w, Du - w), (0,0) 

<a ll'rjl^ + (3 ||cr| 

f a VIV > a \\Du - w\\j^ + (r - {Du - w),C) , Vr G H^fl)*, 

\ 0 IVIlLP(n) — 0 IVllLP(n) + (a-n;,0),VaGHVV* 

f a VIV ^ “ 11^^ “ ^1V + V “ V) , Vr G M{f2), 

[ P IVIlLP(n) ^ P I VllLP(n) + ((^ - 0) , Vct G LP(V 

e ad IHV ~ 

0 e /39|HlLP(a) (V, 

and using (A.3) we can simplify the last expressions with: 

(4.15) ^ G aSgn{Du — w), 


10 


□ 








and 


(4.16) 


0 € {0 e : 


ll«'IILp(ii) 


L^(Q) 


< /?} if w = 0^ 


if w ^ 0. 

Indeed, the norm is an one homogeneous functional and thus its sub differential reads: 

^IMIlp(Q) H ^ {ze (L^’(f2))* : {z,w) = ||n;||LP(Q), {z,(t) < ||cr||Lp(Q), Vcr e L^’(f2)| . 

Clearly, for = 0, the above expression reduces to ||cr||LP(^^) > (^ 5 ^) ? ^ L^(ri) which is 

valid for any z G with ||^||Lg(^^) < 1, i.e., the unit ball of L^(ri). If ^ 0 then the 

subdifferential reduces to the Gateaux derivative of the norm, i.e., d {w) = . 

Finally, from (4.14) we have for every (w^u) G y* 

F|(w,n) > F^iw^u) + {{w,u) - {w,u)) 

^ [ if -u)‘^dx>i j if - uf dx+a^-cp, iw-w,u- u)) 4^ 

i [ if -u)‘^dx>l- [ if - u)‘^ dx + (i^ - (f),w - w) + u-u) 4^ 

Z JQ Z Jq 

f (^-4),w-w} <0,Vtd e iH^ifi))*, 

1 5 fnif - u-ii) <\ J^if - iif, Vn G (L2(11))* 

r 

\ ^d(^l\\f --Wfj iu) ^u-f. 

Combining all the above results, we obtain the optimality conditions (4.11) and (4.12). □ 


Remark 4.5. IFe observe that ifw = 0 then the eonditions (4.11) eoineide with the optimality 
eonditions for the L^-TV minimisation problem (ROF) with parameter a, i.e., 

(4.17) min i\\f-u\\l +a \\Du\\^, 

nGBV(S2) Z 

see also [26]. On the other hand when w ^ 0, the additional eondition (4.12) depends on the 
value of p and as we will see later it allows a eertain degree of smoothness in the final solution 

u. 


4.2. Structure of the solutions. The optimality conditions (4.11) and (4.12) are an 
important asset since we can determine exactly the structure of the solutions for the problem 
{V) as this is determined by the regularising parameters a, (3 and the value of p. 

We initially discuss the cases where the solution u of (V) is a solution of a corresponding 
ROF minimisation problem i.e., = 0. 

Proposition 4.6 (ROF-solutions). Let q be the eonjugate exponent of p ^ (l^cx)] as this is 
defined in (A.2). If 

(4.18) ^>|f2|t, 

a 

then the solution u of {V) eoineides with the solution of the ROF minimisation problem (4.17) 
and w = 0. 


Proof. Let (r(;*,i^*) be a solution pair for (P), then for every {w,u) G L^(fi) x BV(i4), 

\ 11/ - «1lL2(a)+« \\Du* - w*\\^+/3 ||w*||LP(n) ^ ^ 11/ - «llL2(n)+« ll^t^llLP(a) • 


11 








Setting le = 0, we get 

(4.19) ^ 11/ - w*||l2(Q) + « \\Du* - w*\\^ + /3 ||^u’'||Lp(o) < ^ 11/ - «IIl2(Q) + « \\Du\\^ . 

Since w G with p G (1, oc], using the inequality (A.l) we have that 

ll^*ll/W = II^*IIl1(^7) — 1^1 ^ ’ 

and using the condition (4.18) we get 

a < a WDu"^ — + a ||ie*||^ 

<a\\Du*-w*\\j^ + ( 3 \n\'^~" 

(4.20) <a\\Du*-w*\\^ + p\\w*\\^,^^y 
From (4.19) and (4.20) we conclude that for every u G BV(fl), 

^ 11/ - «1Il2(q) + « \\Dn*\\^ ^ ^ 11/ - ^Wlnn) + « > 

i.e., I/* solves 

1 2 1 2 
ug’bv^q) 2 “ “IIl2(Q) + « \\Du\\^ = „g“1?(0) 2 ~ “llL2(n) + « “ M\m + ^ ll^llLP(f2) • 

w=0 

□ 


In fact what we have essentially proved above is that when the condition (4.18) holds then 

TVL^,^(//) = a\\Du\\M. yu G BV(fl). 

Notice also that when (4.18) holds then we can show that le = 0 is an admissible solution but 
in general we cannot prove that this solution is unique. 

The condition (4.18) is valid for any dimension d > 1. It provides a rough threshold for 
obtaining ROF type solutions in terms of the regularising parameters (a,/3 and the image 
domain fl. However, the condition is not sharp since as we will see in the following sections 
we can obtain a sharper estimate for specific data /. 

The following proposition in the spirit of [6, 24] gives more insight into the structure of 
solutions of {V). 

Proposition 4.7. Let f G BV(fl) and suppose that (w^u) G L^(fl) x BV(fl) is a solution 
pair for {V) with p G (1, oc). Suppose that u> f (or u < f ) on an open interval I C Lt then 
(Du — re) [/ = 0 i.e., u' — w on I and \D^u\{I) = 0. 

The above proposition is formulated rigorously via the use of preeise representatives of 
BV functions, see [1], but for the sake of simplicity we rather not get into the details here. 
Instead we refer the reader to [6,24] where the analogue propositions are shown for the TGV 
regularised solutions and whose proofs are similar to the one of Proposition 4.7. 

We now consider the case where the solution is constant in fl, which in fact coincides with 
the mean value / of the data f: 

(4.21) /:= argmin \\\f - u\\l 2 (a) = f f dx. 

u constant ^ do 

Proposition 4.8 (Mean value solution). If the following eonditions hold 

a > 11 / - /||li(Q), 

then the solution of {V) is eonstant and it is equal to /. 



12 


Proof. Clearly, if is a constant solution of (7^), then Du — 0 and from (2.2) we get 
TVL^ ^{u) = 0. Hence, we have u — f. In order to have u — f., from the optimality conditions 
(4.11) and (4.12), it suffices to find a function (j) G HQ(i4) such that (/)(a) = 0(6) = 0 and 

Letting (p{x) = — /) ds, then obviously (/>(a) = 0(6) = 0 and 

|0(a;)| < [ |/(s)-/|c^s < ||/-/||li(Q) < oo- 

J a 

Therefore, ||0||^ < ||/ — / ||li(^)* ^Iso, since L^(i4) C L^(i4) we obtain 

||0|lL9(a) < 1^1 ’ Halloo < \^\h\f - /|Il1(Q)- 

Hence, it suffices to choose a and /3 as in (4.22). □ 



11 / “ fWh^n) 


Fig. 1. Characterisation of solutions of {V) for any data /: The blue/red areas 
correspond to the ROF type solutions {w = 0) and the purple area corresponds to 
the TVL^ solutions (re 7 ^ 0) for 1 < p < oc. We note that the red/purple areas are 
potentially larger/smaller as the conditions we have derived are not sharp. 

In Figure 1, we summarise our results so far. There, we have partitioned the set {a > 
0,/3 > 0 } into different areas that correspond to different types of solutions of the problem 
(V). The brown area, arising from thresholds (4.22) corresponds to the choices of a and /3 that 
produce constant solutions while the blue area corresponds to ROF type solutions, according 
to threshold (4.18). Therefore, we can determine the area where the non-trivial solutions are 
obtained i.e., le 7 ^ 0, see purple region. Note that since the conditions (4.18) and (4.22) are 
not sharp the red and the purple areas are potentially larger or smaller respectively than it 
is shown in Figure 1. 

The following proposition reveals more information about the structure of solutions in the 
case w ^ 0. 

Proposition 4.9 (TVL^-solutions). Let f G BV(i4) and suppose that (w^u) G L^(i4) xBV(i4) 
is a solution pair for {V) with p G (1, oc) and w ^ 0. Suppose that u > f (or u < f) on an 
open interval I C Lt then the solution u of {V) is obtained by 

(4.23) — C{\u'{x)\^~‘^u'(x))' + u{x) = /(x), Vx G / where C — - ^ . 

II^IIlp(q) 


13 












Proof. Since 1 < p < oo^ w ^ 0 using Proposition 4.7 and the second optimality condition of 
(4.12), we have that 


(/) = (3 


Hence, by (4.11) we obtain (4.23) where C 


\u'\P 


\w 


\p-l 

Ilp(O) 


^ 


iip-i • 

IIlp(q) 


□ 


4.3. Exact solutions of {V) for a step function. In what follows we compute explicit 
solutions of the TVL^ denoising model {V) for the case p = 2 for a simple data function. We 
define the step function in L), L > 0 as: 


0 if X G (-iv,0], 

h if X G (0, L). 

We first investigate conditions under which we obtain ROF type solutions, that is rc = 0. 



4 . 3 . 1 . ROF type solutions. We are initially interested in solutions that respect the discon¬ 
tinuity at X = 0 and are piecewise constant. From the optimality conditions (4.11)-(4.12), it 
suffices to find a function v G HQ(fl) such that 


(4.25) ^{-L) = ^{L) = 0, 11011^ < a, m = 


and it is also piecewise affine. It is easy to see that by setting (j){x) = ^(L — |x|), the 
conditions (4.25) are satisfied and the solution u is piecewise constant. The first condition 

of (4.12) implies that ||0||Lg(^) < /3 ^ > (^^)^ and provides a necessary and sufficient 

condition that need to be fulfilled in order for u to be piecewise constant, that is to say 


(4.26) 


u{x) 



if X G (—L, 0] 
if X G (0, L) 


o, \q 1J 


A special case of the ROF-type solution is when u is constant, i.e., when u — f the mean 
value of /. 

We define (^{x) = |(L — |x|) and in that case we have that ||0||^ < a a > ^ and 
||0 |Il9(q) I3> This implies that 


(4.27) 


u — f — — 
■' 2 


a > — and 
“ 2 


/5> n 


h /2L9+i\ ■J 




with —^— = 1. 
p q 


Using now (4.26)-(4.27) we can draw the exact regions in the quadrant of {a >,/3 > 0} that 
correspond to these two types of solutions, see the left graph in Figure 3 for the special case 
p — 2. Notice that in these regions w — 0 and the estimates are valid for any p G (1, oc). 


4 . 3 . 2 . TVL^ type solutions. For simplicity reasons, we examine here only the case p = 2 
with w ^ 0 in Q. However, we refer the reader to Section 6.2 where we compute numerically 
solutions for p ^ 2. Using Proposition 4.9, we observe that the solution is given by the 
following second order differential equation: 

(4.28) — Cu\x) — uix) = /(x), subject to C = ^ -. 

II^IIl2(o) 

Even though we can tell that the solution of (4.28) has an exponential form, the fact that the 
constraint on C depends on the solution rc, creates a difficult computation in order to recover 
u analytically. In order to overcome this obstacle, we consider the one dimensional version of 
the 2-homogeneous analogue of {V) that was introduced in Section 3: 

«eBV(a) ^ ~ ^llL2(a) + « W^u - w\\^ + ||w^||L2(n) • 

k;GL^(0) 
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Similarly to Section 4.1, one can derive the optimality conditions for (4.29). A pair (w^u) is 
a solution of (4.29) if and only if there exists a function (j) G HQ(fl) such that 

(j)' ^u- f, 

(4.30) (j) G aSgn{Du — re), 

0 ^2—hom'^- 


In order to recover analytically the solutions of (V) for p = 2 and determine the purple 
region in Figure 1 it suffices to solve the equivalent model (4.29) where w ^ 0. We may restrict 
our computations only on / = (—A, 0] C and due to symmetry the solution in / = (0, L) is 
given by u{x) + u{—x) = h. The optimality condition (4.30) results to 


(4.31) 


— u'\x) + ku{x) = 0 , where 


1 


and X G / = (—L,0] C fl. 


Then, we get i/(x) = cie^^ + C 2 e“^^ with (/)(x) = + for all x G (—L,0]. Firstly, 

we examine solutions that are continuous which due to symmetry much have the value | at 
the X = 0, i.e., i/(0) = |. Since 0 G Hq(—L,L), we have 0(—T) = 0 and also u\—L) = 0. 
Finally, we require that 0(0) < a. After some computations, we conclude that 


(4.32) 


u{x) 


c^^kx kx if X G (—L,0], tanh(A:L) 2a 

h — cie~^^ — C2e^^ if x G (0, L) k /i ’ 


where ci = 026 ^^^, C 2 = 2 (e^kL^i) &ndk = 

On the other hand, in order to get solutions that preserve the discontinuity at x = 0, we 
require the following: 


(4.33) 

Then we get 

(4.34) u{x) 


0(-L) = 0, u'{-L) = 0, 

?i( 0 ) < ^, 0 ( 0 ) = a. 


cie^^ + C2e 
h — cie~^^ — C2e^^ 


if X G (—L,0], tanh(A:L) 2a 

if X G (0, L) k /i ’ 


where ci = C 2 e^^^, C 2 = ^ 2 kL_i ^tnd k — Notice that the conditions for a and 0 in (4.32) 
and (4.34) are supplementary and thus only these type of solutions can occur, see the quadrant 
of {(a > 0,0 > 0} as it presented in Figure 2. Letting ^(0) = v^f 2 inh(;^), if ^(0) < ^ 
then the solution is of the form (4.32), see the blue region in Figure 2 . On the other hand 
in the complementary green region we obtain the solution (4.34). For extreme cases where 
0 ^ 00 , i.e., /c ^ 0 we obtain ^ which means that there is an asymptote of g 

at (a = ^. Although, we know the form of the inverse function of the hyperbolic tangent, 
we cannot compute analytically the inverse f~^. However, we can obtain an approximation 
using a Taylor expansion which leads to 


(4.35) 





^ hL^ 

h ^ 3(/iL —2(a)’ 


where (a > 0 and a ^ 

Finally, we would like to describe the solution on the limiting case 0 ^ oo. Letting 0 ^ oc 
in (4.32), we have that ci,C 2 ^ | and u{x) | for every x G fl, which in fact is the mean 
value obtained from {V). For the discontinuous solutions, we have that ci,C 2 ^ ^ and 


f - 

u(x) \ ^ 


if X G (-L, 0], 
if X G (0, L), 
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^ 2 —horn 



hL 


Fig. 2. Characterisation of solutions of (4.29) for data / being a step function. The 
green region corresponds to solutions that preserve the discontinuity at x = 0, (4.34), 
while the blue region corresponds to continuous solutions, (4.32), both having an 
exponential form. 


i.e., we converge to the solution (4.26). We also get that 


(4.36) 


w{x) = kC 2 


^2kL-\-kx _ ^—kx 

^2kL—kx _ ^kx 


if X G ( —L, 0], 
if X G (0, L], 


with ||'?^||l2(o) — C 2 k\^e^^{sinh{2kL) — 2kL)^ and C 2 is given either from (4.32) or (4.34). 
Then, in both cases we have w ^ 0 a,s k ^ 0. Observe that the product of P2-hom II'^IIl2(o) 
is bounded as 132-hom cxo for both types of solutions and in fact corresponds to the bounds 
found in (4.26) and (4.27). Indeed, since 


(sinh(2/cL) — 2kL)‘2 



as k —y 0, 


if (a > ^ then 


while if (a < ^ 


^2 —hom ll'*^llL2(r2) 


h 

2 V 


as 1^2—hom ^ OO) 


h2 -hom ||'^|Il2(0) ^ 



as (32—hom ^ CXD- 


The last result is yet another verification of Theorem 3.2 and it shows that there is an 
one to one correspondence, 132-hom II'^IIl2(o) ^ l3i-hom and the purple region of Figure 3 is 
characterised by the solutions obtained in Figure 2. 


5. An image decomposition approach 

In this section, we present another formulation for the problem (P), where we decompose 
an image into a BV part (piecewise constant) and a part that belongs to (smooth). 

Let 1 < p < 00 and 14 C and consider the following minimisation problem: 

^ 11-^ “ “ “ ^IIl2(q) + « \\Du\\^ + 13 || V?;|lLp(f2). 

vew^^p(Q) 

In this way, we can decompose our image into two geometric components. The second term 
captures the piecewise constant structures in the image, whereas the third term captures the 
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^1—hom /^2—/iom| I'^l |l2(0) 


Fig. 3. Characterisation of solutions of (7^) for p = 2 for data / being a step function. 
The type of solutions in the purple region of the left graph are exactly the solutions 
obtained for the 2-honiogenous problem (4.29), on the right graph. 


smoothness that depends on the value of p. In the one dimensional setting, we can prove that 
the problems {V) and (5.1) are equivalent. 

Proposition 5.1. Let = (a, 5) C M, then a pair ('u*,i/*) G W^’^(f2) x BV(ri) is a solution 
of (5.1) if and only if + '?;*) G L^(ri) x BV(ri) is a solution of (V). 

Proof. Let u — u + v then, we have the following 

{v*,u*)e argmin ^ \\f - u - v\\l 2 ^^^ + a\\Du\\^ + (3 ^ 

■ueBV(a) ^ 

^GWI’P(Q) 


{v*,u*)e argmin \ \\f - u - v\\l 2 (^m + a sup {(u, (/>')}+^ ||Vu||Lpm) 
■ueBV(a) ^ 4>ec^(n) 

vew^’P(n) ll-/'lloo<i 

{v*,u*)e argmin L|/- i^||L 2 (Q) « sup {(u, <?i') + (Vu, 0)} +/? || VuIIlp^) 
«gbv(Q) ^ <pec^{n) 

pewi'P(a) \\4>\\oo<^ 

{w*,v*)e argmin I \\f - u\\l 2 ^a) + - w\\j^ + P \\w\\^p,^). 

ueBV(a) ^ ^ ' 

w='Vv 
pGWI’P(Q) 


However, we can eliminate the last constraint since 


(5.2) {w G L^(5^) : 3^ G W^’^(J3), w = Vv}= L^(53). 

Indeed, let w G L^(i2) C L^(i2) for p G (1, oc] and define v{x) = w{s) ds for x G 14 C M. 
Clearly, P — w a.e and by Jensen’s inequality 


\v{x)\^ = 


rx P rx 

/ w{s)ds < / \w{s)\^ ds < oo^ 

J a J a 


and V G W^’^(14) for p G (l,oo). Finally, for the case p = oo, let C > 0 be a constant such 
that \w{x)\ < C a.e. on 14. In that case we have |^’(x)| < |'^(5)| ds < C\Q\ < oo and 

ry 

\v{y) — v{x)\< / \w{s)\ds < C\y — x\^ Vx,2/G14, 

J X 
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i.e., V G from Rademacher’s theorem. Therefore, 

_ 1 _ 2 - 

{w\ u^) G argmin - ||/ - a\\Du-w\\^+ , 

uGBV(O) ^ 
weLP(n) 

where ix* = i/* + x;* and le* = Vv*. □ 

Even though for d = 1 it is true that every function can be written as a gradient, this is 
not true for higher dimensions. In fact, as we show in the following sections, this constraint 
is quite restrictive and for example the staircasing effect cannot always be eliminated in the 
denoising process, see for instance Figure 18. 

The existence of minimisers of (5.1) is shown following again the same techniques as in 
Theorem 2.5. Moreover, due to the strict convexity on the fidelity term of (5.1), one can 
prove that the sum u -\- v ^ BV(i2) is unique for a solution (u^v) G W^’^(i2) x BV(ri). This 
result coincides with the uniqueness of (V) problem for u. Finally, if (xxi,x;i), {u 2 ^V 2 ) are two 
minimisers of (5.1), then from the convexity of T(xx, v) we have for 0 < A < 1 

L{\{ui,Vi) + (1 - \){U 2 ,V 2 )) < \L{ui,Vi) + (1 - \)L{u 2 ,V 2 ). 

Since (xxi,x;i), (xx2,'^2) are both minimisers, the above inequality is in fact an equality. Since 
ui+ vi — U 2 -\- V 2 ^ we obtain 

a\D{\ui + (1 — A)?i 2 |(f^) + P ||V(Ax;i + (1 — ^)'^2 )||lp(o) 

(5.3) = a{\\Durm + (1 - \)\Du2\) + /3(A ||+ (1 - A) ||Vx;2||lp(o)). 

If we assume that 

||V(Ax;i + (1 — A)'^2 )||lp(0) ^ ^ II^'^iIIlp(O) + (1 ~ II^'^2||lp(0) ’ 

then we contradict the equality on (5.3). Hence, the Minkowski inequality becomes an equality 
which is equivalent to the existence of /x > 0 such that Vv 2 = jiVvi. In other words, we have 
proved the following proposition that was also shown in [18] in a similar context: 

Proposition 5.2. Let (xxi, x;i), ( 1 x 25 '^ 2 ) be two minimisers of (5.1). Then 

(5.4) Ui + Vi = U2 + V2, 

(5.5) 3/x > 0 sueh that Vv 2 = fiVvi. 


6 . Numerical Experiments 

In this section we present our numerical simulations for the problem {V). We begin with 
the one dimensional case where we verify numerically the analytical solutions obtained in 
Section 4.3. We also describe the type of structures that are promoted for different values of 
p. Finally, we proceed to the two dimensional case where we focus on image denoising tasks 
and in particular on the elimination of the staircasing effect. 

We start by defining the discretised version of problem {V) 

(6.1) min i|l/-«||2 + TVL^ (^). 

Here TVL^ . : ^ M is defined as 

(6.2) argmin a\\\/u — w\\-^ + f3 \\w\\ ^ 

’ K;G(M’^Xm)2 
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where for x G we set ||x||^ = for x = (xi,X2) G we 

define 


(6.3) 



n^m 


Y1 {ViMiJW + {X2{i,j))'^y 




1 

P 


We denote by V = (Vi, V 2 ) the discretised gradient with forward differences and zero Neu¬ 
mann boundary conditions defined as 


{'Viu)ij — 
i^2u)ij = 


0 

0 


if 1 < i < n, 1 < j < m, 
if z = n, 1 < j < m, 

if 1 < z < n, 1 < j < m, 
if 1 < z < n, j = m. 


where t denotes the step size. The discrete version of the divergence operator is defined as 
the adjoint of V. That is, for every w = (zz;i,zz;2) ^ ^ ^ ^nxm^ have that 

(—divzz;,zz) = (zz;,Vzz) with 


(6.4) 


{diYw)ij 


t 

t 

wi{i,j-l) 
< t 


+ 


t 

t 

W2{i-lJ) 

t 


if 1 ^ j < 1 < 

if j = 1, 1 < z < n, 
if j = m, 1 < z < n, 

if 1 < z < n, 1 < j < m, 
if z = 1, 1 < j < m, 
if z = m, 1 < j < m. 


We solve the minimisation problem (6.1) in two ways. The first one is by using the CVX 
optimisation package with MOSEK solver (interior point methods). This method is efficient 
for small-medium scale optimisation problems and thus it is a suitable choice in order to 
replicate one dimensional solutions. On the other hand, we prefer to solve large scale two 
dimensional versions of (6.1) with the split Bregman method [15] which has been widely used 
for the fast solution of non-smooth minimisation problems. 


6.1. Split Bregman for L^—TVLP. In this section we describe how we adapt the split 
Bregman algorithm to our discrete model (6.1). Letting 2 ; = Vzz — ze, the corresponding 
unconstrained problem becomes 


J_ 2 

(6.5) min - 11/— zzll^-h cr lUII-, -h/3 ||zz;|L , such that z = Vu — w. 

ueR^xm 2 ^ 

^^(I^nxm )2 

^^(I^nxm )2 

Replacing the constraint, using a Lagrange multiplier A, we obtain the following uncon¬ 
strained formulation: 


(6.6) mm - / - u\\l -h a z N -h /3 \\w\\ -h - 2 ; - Vzz + w\\l . 

^^(]^nxm)2 

^^(]^nxm )2 

The Bregman iteration, see [23], that corresponds to the minimisation (6.6) leads to the 
following two step algorithm: 
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(6.7) 

(6.8) 


zk+l 


w;^+^) = argmin ^ ||/ - u\\l + a ||z||;^ + P \\w\\p + ^ 

u,z,w ^ ^ 

j^k+1 ^Ijk _ y^fc+1 _ ^k+i^ 


b'^ - z + Vu - w 


5 

2 


Since solving (6.7) at once is a difficult task, we employ a splitting technique and minimise 
alternatingly for z and w. This yields the split Bregman iteration for our method: 


(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

( 6 . 12 ) 




= argmin - 11/- «||2 + - 

Z Z 

. II II ^ 

z ^ = argmin a \\z\\-^ + — 

zG(M’^Xm)2 2 

. ^11 II ^ 

= argmin p||ic||^+ — 

K;G(M’^Xm)2 2 

J^k+l ^ + ^fc+l _ y^fc+l _ y^/o+1 


+ z- + 




^k ^ ^fc+l _ y^fc+l ^ 


2 

2 

5 

2 

W 


2 

5 

2 


Next, we discuss how we solve each of the subproblems (6.9)-(6.11). The first-order opti¬ 
mality condition of (6.9) results into the following linear system: 

(6.13) {I - XA)u = f - Adiv(6'= + z’^ - w^) . 

^ -V-^ -V-^ 

A C 

Here A is a sparse, symmetric, positive definite and strictly diagonal dominant matrix, thus we 
can easily solve (6.13) with an iterative solver such as conjugate gradients or Gauss-Seidel. 
However, due to the zero Neumann boundary conditions, the matrix A can be efficiently 
diagonalised by the two dimensional discrete cosine transform. 


(6.14) 


A = W^^DWn 


where here Wnm is the discrete cosine matrix and D = diag^jii^ • • • is the diagonal 

matrix of the eigenvalues of A. In that case, A has a particular structure of a block symmetric 
Toeplitz-plus-Hankel matrix with Toeplitz-plus-Hankel blocks and one can obtain the solution 
of (6.9) by three operations involving the two dimensional discrete cosine transform [16] as 
follows: Firstly, we calculate the eigenvalues of A by multiplying (6.14) with ei = (1,0, • • • , 0)"^ 
from both sides and using the fact that WnmWnm — WnmWnm = Inm-> we get 


(6.15) 


\WnYYlAC\\i 

\Wnm^\\i 


i = 1, 2, • • • , nm. 


Then, the solution of (6.9) is computed exactly by 


(6.16) 


u = WJ^D-^Wr, 


c. 


The solution of the subproblem (6.10) is obtained in a closed form via the following shrinkage 
operator^ see also [15,30]. Indeed, for i = 1,2 we have 


(6.17) 


= shrinka (6f - + wf) 


:= max 


9i 


9i 


Finally, we discuss the solution of the subproblem (6.11). In the spirit of [29], we solve (6.11) 
by a fixed point iteration scheme. Letting = y and g = —b^ — z^^^ + the first-order 

optimality condition of (6.11) becomes 


(6.18) 


^ + w -g = 0 
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For given we obtain by the following fixed point iteration 

II hWP—^ 

Vi \\w \ 

(6.19) 


= 


K.\W 


k\p-2 \\w^\\P~^ ’ 


under the convention that 0/0 = 0. We can also consider solving the p-homogenous analogue 
{Vp-hom)^ where for certain values of p, e.g. p = 2, we can solve exactly (6.19), since in that 


,,^ + 1 _ Pi 
i K-\-\ 

difference between these two methods. 


case W- ' ^ However, we observe numerically that there is no significant computational 


Let us finally mention that since we do not solve 
exactly all the subproblems (6.9)-(6.11), we do not have a convergence proof for the split 
Bregman iteration. However in practice, the algorithm converges to the right solutions after 
comparing them with the corresponding solutions obtained with the CVX package. 


6.2. One dimensional results. For this section, we set m = 1 and thus u G 
w G Initially, we compare our numerical solutions with the analytical ones, obtained 

in Section 4.3 for the step function, setting p — 2^ h — 100, L — 1 and = [—1,1]. The 
domain is discretised into 2000 points. We first examine the cases of where ROF solutions 
are obtained, i.e., the parameters a and (3 are selected according to the conditions (4.26) and 
(4.27), see Figure 4. There we see that the analytical solutions coincide with the numerical 
ones. 






— TVL^ numerical 
» TVL^ exact 

— step function 






(c) (ROF) : a = 60, /? = 1300 


Fig. 4. Comparison between numerical solutions of {V) and the corresponding ana¬ 
lytical solutions obtained in Section 4.3. The parameters a and P are chosen so that 
conditions (4.26) and (4.27) are satisfied. 


Now, we proceed by computing the non-ROF solutions. The numerical solutions are solved 
using the 2-homogeneous analogue of (4.29), since we have proved that the 1-homogeneous and 
p-homogeneous problems are equivalent modulo an appropriate rescaling of the parameter /3, 
see Proposition 3.2. In fact, as it is described in Figure 3, in order to obtain solutions from the 
purple region, it suffices to seek solutions for the 2-homogeneous (4.29). Notice also that these 
solutions are exactly the solutions obtained solving a Huber TV problem, see Proposition 3.3. 
The analytical solutions are given in (4.32) and (4.34) and are compared with the numerical 
ones in Figure 5, where we observe that they coincide. We also verify the equivalence between 
the 1-homogeneous and 2-homogeneous problems where a is fixed and [3 is obtained from 
Proposition 3.2, see Figure 5c. 

We continue our experiments for general values of p focusing on the geometric behaviour 
of the solutions as p increases. In order to compare the solutions for p G (l,oo), we fix 
the parameter a and choose appropriate values of /3 and p. We choose a and /3 so that they 

belong to the purple region in Figure 3, i.e., [3 < and [3 < |( ^^^^^ ) g, hence non-ROF 

solutions are obtained. We set p = {|, |,2,3,4,10} and for the solutions that preserve the 
discontinuity we select (3 = {72,140,430,1350,2400,6800} with fixed a = 20 (observe that 

P < (^^)^<^ is valid in any case), see Figure 6a. For the continuous cases, we set a = 60 
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(a) TVL^ : a = 20, 

f^2 — hom ~ 450 


(b) TVL^ : Of = 60, 
^2 —horn ~ 450 


(c) 1 and 2-homogeneous: 
(X = 15, /52 —/iom~450, 

^1 — hom ~ ^2 — hom ||'^||2 


Fig. 5. Comparison between numerical and analytical solutions obtained in Section 
4.3 for the step function, by solving the 2-homogeneous problem (4.29). The parame¬ 
ters a and (3 are chosen so that conditions (4.32) and (4.34) are satisfied. The last plot 
indicates the equivalence between the 2-homogeneous (4.29) and the 1-homogeneous 
problem {V). 


and (3 — {50,110,430,1700,3000,9500} (again the conditions a > ^ and /3 < ^ 

hold), see Figure 6b. We observe that for p = the solution has a similar behaviour to 
p = 2^ but with a steeper gradient at the discontinuity point. Moreover, the solution becomes 
almost constant near the boundary of f]. On the other hand, as we increase the slope of 
the solution near the discontinuity point reduces and it becomes almost linear with a relative 
small constant part near the boundary. 



Fig. 6. Step function: The types of solutions for the problem (V) for different values of p. 



Fig. 7. Piecewise afhne data: TVL^^ solution with o = 1, /3 = 620. 
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The linear structure of the solutions that appears for large p motivates us to examine the 
case of a piecewise linear data / defined as 


( 6 . 20 ) 


fix) 


Xx if X G (—T, 0], 

Ax + /i if X G (0, L], 


see Figure 7. We set again 12 = [—1,1], A = ^ and the data are discretised in 2000 points. 
As we observe, the reconstruction for p = 15 behaves almost linearly everywhere in 12 except 
near the boundary. In the follow up paper [7], where the case p = oc is examined in detail, 
the occurrence of this linear structure is justified. 

In the last part of this section, we discuss the image decomposition approach presented 
in Section 5. We treat a more complicated one dimensional noiseless signal with piecewise 
constant, affine and quadratic components and solve the discretised version of (5.1) using 
CVX under MOSEK. We verify numerically the equivalence between (5.1) and {V) for p = 2, 
i.e., (Vx,i/ + x) corresponds to (w^u) where (v^u) and {w^u) are the solutions of (5.1) and {V) 
respectively, see Figure 8. We also compare the decomposed parts u^v for two different values 
of p (| and 10). In order to have a reasonable comparison on the corresponding solutions, 
the parameters (a,/3 are selected such that the residual \\f — u — vW^ is the same for both 
values of p. As we observe, the v decomposition with p = | promotes some flatness on the 
solution compared to p = 2, compare Figures 8b and 9a. On the other hand for p = 10, the 
V component promotes again almost affine structures. Figure 9b. Notice, that in both cases 
the V components are continuous. In fact, this is confirmed analytically for every 1 < p < oo, 
since in dimension one W^’^(i2) C (7(12). 



Fig. 8. Numerical results on the image decomposition approach (5.1) for p = 2, see 
Section 5. 


6.3. Two dimensional results. In this section we consider the two dimensional case where 
u G w G (M’^x’^)2 with m > 1 and 12 denotes a rectangular/square image domain. We 

focus on image denoising tasks and on eliminating the staircasing effect for different values of 
p. We use here the split Bregman algorithm proposed in Section 6.1. 
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(a) Decomposition of the data in Fig- (b) Decomposition of the data in Fig¬ 
ure 8a for p = I ure 8a for p = 10 


Fig. 9. Decomposition of the data in Figure 8a into u^v parts for p = | and p = 10. 

The value P = | produces a v component with flat structures while p = 10 produces a 
component with almost affine structures. In both cases we have \\f — u — vW^ = 6.667. 

We start with the image in Figure 10, i.e., a square with piecewise affine structures. The 
image size is 200 x 200 pixels at a [0,1] intensity range. The noisy image, Figure 10b, is a 
corrupted version of the original image. Figure 10a, with Gaussian noise of zero mean and 
variance a — 0.01. 



(a) Square (b) Noisy square: PSNR=20.66 and 

SSIM=0.1791 

Fig. 10. Square with piecewise affine structures and its noisy version with a = 0.01. 



In Figure 11, we present the best reconstructions results in terms of two quality measures, 
the Peak Signal to Noise Ratio (PSNR) and the Struetural Similarity Index (SSIM), see [31] 
for the definition of the latter. In each case, the values of a and /3 are selected appropriately 
for optimal PSNR and SSIM. Our stopping criterion is the relative residual error becoming 
less than 10“^ i.e.. 


( 6 . 21 ) 


\U 


k+1 


— U 


/ fc+11 


< 10 


-6 


Finally, for computational efficiency, we fix A = lOa when 1 < p < 4 and A = lOOOa when 
p > 4 (empirical rule). We observe that the best reconstructions in terms of the PSNR have 
no visual difference among p = |, 2 and 3 and staircasing is present. Figures 11a, 11b and lie. 
This is one more indication that the PSNR - which is based on the squares of the difference 
between the ground truth and the reconstruction - does not correspond to the optimal visual 
results. However, the best reconstructions in terms of SSIM are visually better. They exhibit 
significantly reduced staircasing for p = | and p = 3 and is essentially absent in the case of 
p = 2, see Figures lid, lie and Ilf. 

We can also get a total staircasing elimination by setting higher values for the parameters 
a and /3, as we show in Figure 12. There, one observes that on one hand as we increase p. 
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(a) TVL5: a = 0.1,/3 = 2.5, (b) TVL^: a = 0.1,/3 = 13.5, (c) TVL^: a = 0.1,/? = 76, 

PSNR=33.63 PSNR=33.68 PSNR=33.70 



(d) TVL5: a = 0.3,/3 = 7.7, (e) TVL^: a = 0.3,/3 = 34, (f) TVL^: a = 0.3,/3 = 182, 

SSIM=0.9669 SSIM=0.9706 SSIM=0.9709 


Fig. 11. Best reconstructions in terms of PSNR and SSIM for p = 2,3. 


almost affine structures are promoted - see the middle row profiles in Figure 12 - and on the 
other hand these choices of a, [3 produce a serious loss of contrast that however can be easily 
treated via the Bregman iteration. 

Contrast enhancement via Bregman iteration was introduced in [23], see also [3] for an 
application to higher-order models. It involves solving a modified version of the minimisation 
problem. Setting vP — /, for A; = 1, 2,..., we solve 

2 

^ -h (a ||V?i - w\\-^+ f3\\w\\p, 


( 6 . 22 ) 


A*-1-1 

u = argmm 

U^]^nxm 

nXm\2 


f + v’^ 


u 


wG 


Instead of solving (6.1) once for fixed a and /3, we solve a sequence of similar problems 
adding back a noisy residual in each iteration which results to a contrast improvement. For 
stopping criteria regarding the Bregman iteration we refer to [23]. In Figure 13 we present 
our best Bregmanised results in terms of SSIM. There, we notice that Bregman iteration leads 
to a significant contrast improvement, in comparison to the results of Figure 12. In fact, we 
observe that the Bregmanised TVL^ (first-order regularisation), can achieve reconstructions 
that are visually close to the second-order Bregmanised TGV^, compare Figures 13e and 
13f. The second-order TGV^ and Bregmanised TGV^ are solved using the Chambolle-Pock 
primal-dual method [10]. 

We continue our experimental analysis with a radially symmetric image, see Figure 14. In 
Figure 15, we demonstrate that we can achieve staircasing-free reconstructions for p = 2, 3 
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and 7. In fact, as we increase we obtain results that preserve the spike in the centre of 
the circle, see Figure 15d. This provides us with another motivation to examine the p — oo 
case in [7]. The loss of contrast can be again treated using the Bregman iteration (6.22). The 
best results of the latter in terms of SSIM are presented in Figure 16, for p = 2,4 and 7 and 
they are also compared with the corresponding Bregmanised TGV^. We observe that we can 
obtain reconstructions that are visually close to the TGV^ ones and in fact notice that for 
p = 7, the spike on the centre of the circle is better reconstructed compared to TGV^, see 
also the surface plots in Figure 17. 

We conclude with numerical results for the image decomposition approach of Section 5 
which we solve again using the split Bregman algorithm. Recall that in dimension two, the 
solutions of (5.1) will not necessarily be the same with the ones of (V). In fact, we observe 
that (5.1) cannot always eliminate the staircasing, see for instance Figure 18. Even though, 
we can easily eliminate the staircasing both in the square and in the circle by applying TVL^ 
regularisation. Figures 18b and 18d, we cannot obtain equally satisfactory results by solving 
(5.1). While using the latter we can get rid of the staircasing in the circle. Figure 18c, this is 
not possible for the square. Figure 18a, where we observe - after extensive experimentation - 
that no values of a and /3 lead to a staircasing elimination. This is analogous to the difference 
between TGV^ and the TV-TV^ infimal convolution of Chambolle-Lions [9]. 

However, as we mentioned before, the strength of the formulation (5.1) lies on its ability 
to efficiently decompose an image into piecewise constant and smooth parts. We show that 
in Figure 19, for the image in Figure 18c. 

7. Conclusion 

We have introduced a novel first-order, one-homogeneous TV-L^ infimal convolution type 
functional for variational image regularisation. The TVL^ functional constitutes a very general 
class of regularisation functionals exhibiting diverse smoothing properties for different choices 
of p. In the case p = 2 the well-known Huber TV regulariser is recovered. 

We studied the corresponding one dimensional denoising problem focusing on the structure 
of its solutions. We computed exact solutions of this problem for the case p = 2 for simple one 




(a) TVL5: a = 1, 

/? = 25, SSIM=0.9391 


(b) TVL^; a = l, 


p = 116, SSIM=0.9433 


(c) TVL^; a = 1, 
p = 438, SSIM=0.9430 


(d) TVL^; a = 2, 

P = 5000, SSIM=0.9001 


Fig. 12. Staircasing elimination for p = |,2,3 and 7. High values of p promotes 
almost affine structures. 
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(a) TV: a = 0.2, SSIM=0.9387 (b) TVL^: a = 1, /3 = 116, (c) TGV^: a = 0.12, /? = 0.55, 

SSIM=0.9433 SSIM=0.9861 



(d) Bregmanised TV: a = 1, (e) Bregmanised TVL^: o; = 2, (f) Bregmanised TGV^: a = 2, 

SSIM=0.9401, 4th iteration j] = 220, SSIM=0.9778, 4th (3 = 10, SSIM=0.9889, 8th 

iteration iteration 


Fig. 13. First Row: Best reconstruction in terms of SSIM for TV, TVL^ and TGV^. 
Second Row: Best reconstruction in terms of SSIM for Bregmanised TV, TVL^ and 
TGV^ 




(a) Gircle 


(b) Noisy circle: SSIM=0.2457 


Fig. 14. Image with symmetric radial structures and its noisy version with a = 0.01. 


dimensional data. Hence, as an additional novelty in our paper we presented exact solutions 
of the one dimensional Huber TV denoising problem. 

Numerical experiments for several values of p indicate that our model leads to an elimination 
of the staircasing effect. We show that we can further enhance our results by increasing the 
contrast via a Bregman iteration scheme and thus obtaining results of similar quality to those 
of TGV^. Furthermore, as p increases the structure of the solutions changes from piecewise 
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(a) TVL5: a = 0.8, 
/? = 17, SSIM=0.8909 


(b) TVL^: a = 0.8, 
/? = 79, SSIM=0.8998 


(c) TVL^: a = 0.8, 

^ = 405, SSIM=0.9019 


(d) TVL^: a = 0.8, 

/? = 3700, SSIM=0.9024 


Fig. 15. Better preservation of spike-like structures for large values of p. 


smooth to piecewise linear and the model, in contrast to TGV^, is capable of preserving sharp 
spikes in the reconstruction. This observation motivates a more detailed study of the TVL^ 
functionals for large p and in particular for the case p = oo. 

This concludes the first part of the study of the TV-L^ model for p < oc. The second 
part [7], is devoted to the p — oo case. There we explore further, both in an analytical and 
an experimental level, the capability of the TVL^ model to promote affine and spike-like 
structures in the reconstructed image and we discuss several applications. 
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(a) Original 


(b) Bregmanised TVL^ (c) Bregmanised TVL^ 


(d) Bregmanised TGV 









(e) Original: 
central part zoom 


(f) Bregmanised TVL^: 
central part zoom 


(g) Bregmanised TVL^: 
central part zoom 


(h) Bregmanised TGV: 
central part zoom 


Fig. 17. Surface plots of the images in Figure 16. Notice how high values of p, here 
for p = 7, can preserve the sharp spike in the middle of the image. 
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P = 120, SSIM=0.9268 p = 70, SSIM=0.8994 


Fig. 18. Gomparison between the model (5.1) for p = 2 and TVL^: Staircasing 
cannot be always eliminated. 
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(a) Piecewise constant (b) Smooth component v (c) Middle row profiles 
component u 


Fig. 19. Geometric decomposition of the image in Figure 18d into a piecewise con¬ 
stant and smooth component, u and v respectively, by solving (5.1). 
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Appendix A. Radon Measures and functions of bounded variation 

In what follows C is an open, bounded set with Lipschitz boundary whose Lebesgue 
measure is denoted by \Q\. We denote by Ad(fl,M^) (and Al(ri) if d=l) the space of finite 
Radon measures on Q. The total variation measure of /i G AI(ri,M^) is denoted by |/i|, 
while we denote the polar decomposition of /i by /x = sgn(/x)|/i|, where sgn(/i) = 1 |/i|-almost 
everywhere. 

Recall that the Radon norm of a M^-valued distribution T on is defined as 

||r||>, sup {{u, cf>) : cf> e ||</>|loo < l} • 

It can be shown that ||T||ai < if and only if T can be represented by a measure fi G 
Ad(ri,M^) and in that case WbWm — ImK^)* 

A function u G L^(fl) is a function of bounded variation if its distributional derivative Du 
is representable by a finite Radon measure. We denote by BV(fl), the space of functions of 
bounded variation which is a Banach space under the norm 

II^IIbv(o) II^IIli(o) + • 

The term 

\\Du\\j^^ sup Mdiv</)da; : 0 € ||<?!>||oo < 1 

is called the total variation of also commonly denoted by T\{u). From the Radon-Nikodym 
theorem, the measure fi can be decomposed into an absolutely continuous and a singular part 
with respect to the Lebesgue measure that is Du = D^u + D^u. Here, D^u — VuC^, i.e., 
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Vu denotes the Radon-Nikodym derivative of D^u with respect to C^. When d — Vu is 
simply denoted by u'. 

We will also use the following basic inequality regarding inclusions of spaces 


(A.l) II^IIlpi(O) — 1 ^ 1 ^^ II^IIl?'2(0) ’ ^ ^ Pi < P2 ^ OO. 

Unless otherwise stated q denotes the Holder conjugate of the exponent p, i.e., 

^ ifpG(l,oc), 
if p = oo. 


(A.2) 



Regarding the sub differential of the Radon norm we have that it can be characterised, at least 
for (7o functions, as follows [6] 


(A.3) d IHI^ (p) n = Sgn(p) n 

where here Sgn(p) denotes the set-valued sign 

(A.4) Sgn(p) = {v ^ L^{n) n L^(54, |p|) : ||^||^ <l,v = sgn(p), |p| - a.e.} . 
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